library(tidyverse)
library(viridis)
library(plotly)

knitr::opts_chunk$set(
    echo = TRUE,
    warning = FALSE,
    fig.width = 8, 
  fig.height = 6,
  out.width = "90%"
)

options(
  ggplot2.continuous.colour = "viridis",
  ggplot2.continuous.fill = "viridis"
)
scale_colour_discrete = scale_colour_viridis_d
scale_fill_discrete = scale_fill_viridis_d
theme_set(
  theme_minimal() + 
    theme(
      legend.position = "bottom",
      plot.title = element_text(hjust = 0.5)
    )  
)

Data Source

IPUMS Health Surveys: NHIS is a harmonized set of data covering more than 50 years (1963-present) of the National Health Interview Survey (NHIS). The NHIS is the principal source of information on the health of the U.S. population, covering such topics as general health status, the distribution of acute and chronic illness, functional limitations, access to and use of medical services, insurance coverage, and health behaviors. On average, the survey covers 100,000 persons in 45,000 households each year. The IPUMS NHIS facilitates cross-time comparisons of these invaluable survey data by coding variables identically across time. Our analysis will use data from 2015 to 2021, which covers the COVID-19 period.

Data cleaning

We pulled out data from IPUMS Health Surveys: NHIS and will limit our analysis using data from 2015 to 2021. To analyze the trend of anxiety prevalence, frequency and level from 2015 to 2021, we will focus on anxiety indicators listed below:

  • WORFREQ:How often feel worried, nervous, or anxious
  • WORRX: Take medication for worried, nervous, or anxious feeings
  • WORFEELEVL: Level of worried, nervous, or anxious feelings, last time

Core demographic and Social economic status indicators listed below are also included in this analysis: - AGE:Age - SEX:Biological sex - MARST:Current marital status - POVERTY:Ratio of family income to poverty threshold

Responses indicate Unknown or not applied are excluded from our analysis.

anxiety = 
  read_csv("data/nhis_data01.csv") %>% 
  janitor::clean_names() %>% 
  filter(year>=2015) %>% 
  select(year, worrx, worfreq, worfeelevl, age, sex, marst, poverty) %>% 
  mutate(
    sex = recode_factor(sex, 
                        "1" = "Male", 
                        "2" = "Female"),
    marst = recode_factor(marst, 
                        "10" = "Married", "11" = "Married", "12" = "Married", "13" = "Married",
                        "20" = "Widowed",
                        "30" = "Divorced",
                        "40" = "Separated",
                        "50" = "Never married"),
    poverty = recode_factor(poverty, 
                        "11" = "Less than 1.0", "12" = "Less than 1.0", 
                        "13" = "Less than 1.0", "14" = "Less than 1.0",
                        "21" = "1.0-2.0", "22" = "1.0-2.0", 
                        "23" = "1.0-2.0", "24" = "1.0-2.0", 
                        "25" = "1.0-2.0",
                        "31" = "2.0 and above","32" = "2.0 and above",
                        "33" = "2.0 and above","34" = "2.0 and above",
                        "35" = "2.0 and above","36" = "2.0 and above",
                        "37" = "2.0 and above","38" = "2.0 and above"),
    worrx = recode_factor(worrx,
                          '1' = "no", 
                          '2' = "yes"),
    worfreq = recode_factor(worfreq, 
                            '1' = "Daily", 
                            '2' = "Weekly", 
                            '3' = "Monthly", 
                            '4' = "A few times a year", 
                            '5' = "Never"),
    worfeelevl = recode_factor(worfeelevl, 
                               '1' = "A lot", 
                               '3' = "Somewhere between a little and a lot", 
                               '2' = "A little")
    ) %>% 
  filter_at(vars(worrx, worfreq, worfeelevl), any_vars(!is.na(.)))

Percentage of people reported taken medication for worried, nervous, or anxious feelings

According to the figure, from 2015 to 2021, the percentage of people who report taking medication for worry, stress or anxiety is constantly increasing from 9.13% in 2015 to 13.57% in 2021. We can observe a rapid increase from 2017 to 2019 and, contrary to our expectations, a relatively slow increase from 2019 to 2020. The effect of COVID-19 on anxiety is not evident in this plot.

anxiety %>%
  drop_na(worrx) %>% 
  group_by(year, worrx) %>% 
  summarize(wor_num = n()) %>% 
  pivot_wider(
    names_from = worrx,
    values_from = wor_num
  ) %>% 
  mutate(
    wor_percentage = yes/(no + yes)*100,
    text_label = str_c(yes, " out of ", no + yes)
  ) %>% 
  ungroup() %>% 
  plot_ly(
    y = ~wor_percentage,
    x = ~year,
    color = ~year,
    type = "bar", 
    colors = "viridis",
    text = ~text_label
  ) %>% 
  layout(
    xaxis = list (title = ""),
    yaxis = list (title = "Percentage"),
    showlegend = FALSE
  ) %>% 
  hide_colorbar()

Stratify by Biological sex, family poverty status and martial status

Biological sex, family income, and marital status have long been recognized as important factors associated with mental health. In the following sections, we will stratify the percentages according to these factors to investigate trends among people with different characteristics.

Stratify the reported percentage of people taking medication for worried, nervous, or anxious feelings by biological sex, we can observe a faster increase in the percentage among females from 14.41% in 2018 to 16.52% in 2019. This increase may indicate COVID-19-induced anxiety in females. Among males, the percentage is relatively stable from 2018 to 2019, while there is an increase from 2017 to 2018 and 2020 to 2021.

anxiety %>%
  drop_na(sex, worrx) %>% 
  group_by(sex, year, worrx) %>% 
  summarize(wor_num = n()) %>% 
  pivot_wider(
    names_from = worrx,
    values_from = wor_num
  ) %>% 
  mutate(
    wor_percentage = yes/(no + yes)*100,
    text_label = str_c(yes, " out of ", no + yes)
  ) %>% 
  ungroup() %>% 
  plot_ly(
    y = ~wor_percentage,
    x = ~year,
    color = ~sex,
    type = "bar", 
    colors = "viridis",
    text = ~text_label
  ) %>% 
  add_trace(
    x = ~year,
    y = ~wor_percentage,
    color = ~sex,
    type='scatter',
    mode='lines+markers'
  ) %>% 
  layout(
    xaxis = list (title = ""),
    yaxis = list (title = "Percentage"),
    legend = list(orientation = 'h')
  )

Stratify by ratio of family income to poverty threshold

Stratify the percentage of people reported taken medication for worried, nervous, or anxious feelings by the ratio of household income to the poverty line, we can see that the percentage among the poorest stratum decreased rapidly from 17.30% in 2017 to 15.84% in 2018, which is the opposite of what happened in the other two strata.Although the percentage of the poorest strata decreased rapidly from 2017 to 2018, they still had the highest percentage of the three strata, and this decrease was followed by a rapid increase from 15.84% in 2018 to 18.58% in 2019, which may indicate that people belonging to the poorest stratum are more susceptible to anxiety caused by COVID-19. From 2020 to 2021, the percentage decreases for the other two strata, while for the richest strata, the percentage steadily increases.

anxiety %>%
  drop_na(poverty, worrx) %>% 
  group_by(poverty, year, worrx) %>% 
  summarize(wor_num = n()) %>% 
  pivot_wider(
    names_from = worrx,
    values_from = wor_num
  ) %>% 
  mutate(
    wor_percentage = yes/(no + yes)*100,
    text_label = str_c(yes, " out of ", no + yes)
  ) %>% 
  ungroup() %>% 
  plot_ly(
    y = ~wor_percentage,
    x = ~year,
    color = ~poverty,
    type = "scatter", 
    mode = "lines+markers",
    colors = "viridis",
    text = ~text_label
  ) %>% 
  layout(
    xaxis = list (title = ""),
    yaxis = list (title = "Percentage"),
    legend = list(orientation = 'h')
  )

Stratify by current martial status

Stratify the percentage of people reported taken medication for worried, nervous, or anxious feelings by current martial status, we can observe a rapid increase from 10.61% in 2017 to 15.60% in 2019 among those widowed, while it is difficult to tell whether this increase is caused by COVID-19 as it starts at 2017. There is also a rapid increase among those who are separated, from 14.31% in 2019 to 17.49% in 2020. Considering the timing, this could be an effect of COVID-19.

anxiety %>%
  drop_na(marst, worrx) %>% 
  group_by(marst, year, worrx) %>% 
  summarize(wor_num = n()) %>% 
  pivot_wider(
    names_from = worrx,
    values_from = wor_num
  ) %>% 
  mutate(
    wor_percentage = yes/(no + yes)*100,
    text_label = str_c(yes, " out of ", no + yes)
  ) %>% 
  ungroup() %>% 
  plot_ly(
    y = ~wor_percentage,
    x = ~year,
    color = ~marst,
    type = "scatter",
    mode='lines+markers',
    colors = "viridis",
    text = ~text_label
  ) %>% 
  layout(
    xaxis = list (title = ""),
    yaxis = list (title = "Percentage"),
    legend = list(orientation = 'h')
  )

Frequency of worried, nervous, or anxious feelings

From this bar plot about how often people feel worried, nervous, or anxious, we can observe that the frequency is steadily increasing from 2015 to 2021. There is also a rapid increase from 2019 to 2020, which could be caused by COVID-19.

anxiety %>% 
  drop_na(worfreq) %>% 
  group_by(year, worfreq) %>% 
  summarize(count = n()) %>% 
  group_by(year) %>% 
  summarize(
     percentage=100 * count/sum(count),
     sum_count = sum(count),
     worfreq = worfreq,
     count=count
  ) %>% 
  mutate(
    text_label = str_c(count, " out of ", sum_count)
  ) %>% 
  plot_ly(
    y = ~percentage,
    x = ~year,
    color = ~worfreq,
    type = "bar", 
    colors = "viridis",
    text = ~text_label
  ) %>% 
  layout(
    xaxis = list (title = ""),
    yaxis = list (title = "Percentage"), 
    barmode = 'stack',
    legend = list(orientation = 'h')
  )

Level of worried, nervous, or anxious feelings

From this bar plot about the level of worried, nervous, or anxious feelings people felt last time, we can observe a relatively large increase from 2018 to 2019 in the number of people who felt worried, stressed, or anxious a lot or between a little and a lot, which could be an effect of COVID-19.

anxiety %>%
  drop_na(worfeelevl) %>% 
  group_by(year, worfeelevl) %>% 
  summarize(count = n()) %>% 
  group_by(year) %>% 
  summarize(
     percentage=100 * count/sum(count),
     sum_count = sum(count),
     worfeelevl = worfeelevl,
     count=count
  ) %>% 
  mutate(
    text_label = str_c(count, " out of ", sum_count)
  ) %>% 
  plot_ly(
    y = ~percentage,
    x = ~year,
    color = ~worfeelevl,
    type = "bar", 
    colors = "viridis",
    text = ~text_label
  ) %>% 
  layout(
    xaxis = list (title = ""),
    yaxis = list (title = "Percentage"), 
    barmode = 'stack',
    legend = list(orientation = 'h')
  )